Je donne une liste des problèmes posés par les étudiants pendant le TP (groupe 2 de 4M074 2018). Comme il s’agit beaucoup de détails pour faire des codes jolis, je trouve que c’est plus raisonable d’apprendre les détails petit à petit. Je sélectionne les problèmes qui sont (sans doute) intéressants pour qu’on peut avancer ensemble.

N’hésitez pas de m’envoyez un email pour poser des questions ou corriger mon français.

misc:

1 TP4(corrigé) et TP5(sujet) 20-02-2018

2 Environnements de R

On travaille principalement sur Jupyter Notebook avec IRkernel. Vous pouvez aussi utiliser Rstudio. Voici un guide pour l'installation de Jupyter et IRkernel pour tous les OS.

Installation de Python et Jupyter

Si vous avez des problèmes sur la configuration, n'hésitez pas de me parler. De plus, sagecloud nous permet de compiler les notebooks en ligne (graduit).

2.1 Conseils pour les débudants de R

Si vous n'avez jamais suivi un cours de R, vous trouvrez dans les liens suivants deux documents par madame Tabea Rebafka pour le cours 4M015.

3 Rappels

3.1 Matrices des figures

Cas simple :

par(mfrow=c(1,3))
hist(rnorm(500), breaks=30)
hist(rexp(500), breaks=30)
hist(runif(500), breaks=30)
Figure 1: Matrice d’histogrammes

Cas (un peu) compliqué:

# One figure in row 1 and two figures in row 2
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
hist(rnorm(500), breaks=30)
hist(rexp(500), breaks=30)
hist(runif(500), breaks=30)
Figure 2: Matrice d’histogrammes

3.2 lines et curve

hist(rexp(300), freq=FALSE, breaks = 40)
x<-seq(0,8,by=.1)
# ici, x est un vecteur
lines(x, dexp(x), col = 'darkred')
Figure 3: ex. de la fonction lines 

hist(rexp(300), freq=FALSE, breaks = 40)
# ici, x est qqch abstrait
# pour aller plus loin  ?expression
curve(dexp(x), col = 'darkred', add=TRUE, from=.001)
Figure 4: ex. de la fonction curve

3.3 list to data.frame

a <- list()
b <- list() 
a$X <- 1
a$Y <- 2
b$X <- 3
b$Y <- 4
data.frame(do.call(rbind, list(a,b)))
##   X Y
## 1 1 2
## 2 3 4

3.4 Visualisation 3D par plotly et ggplot2

library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(ggplot2)

a <- seq(-10,10,0.1)
b <- a

f <- function(x,y){
      return(exp(-.01 *((x-1)^2 + (x+1)^2 + y^2)))
}

# z est une matrice t.q. z(i,j) = f(a_i,b_j) où
# (a_i,b_j)est le maillage. 

### par plotly
plot_ly(x = a, y = b, z = sapply(1:length(b), function (i) return(f(a,b[i]))), type = 'surface')

ggplot2_plotly&Figure 5: Visualisation par ggplot2 

### par ggplot2
A = length(a)
B = length(b)
df <- data.frame(a=rep(a,times=B),b=rep(b,each=A),f_z=c(sapply(1:length(b), function (i) return(f(a,b[i])))))
ggplot(df, aes(a, b, z=f_z)) + geom_contour(bins = 60, aes(colour = ..level..))
Figure 5: Visualisation par ggplot2

4 Remarques du TP

4.1 TP1

4.1.1 un exemple de optimize

On remarque que cette fonction n’est pas (du tout) magique pour l’optimisation non-convexe. Mais dans le cas convexe, il marche très bien.

f<-function(x){return(x^2+sin(5*x))}
opt<-optimize(f, interval=c(-3,3), maximum=F, tol=0.01)
x_min<-opt$objective
curve(f(x), from=-3, to=3)
points(x_min, f(x_min), col='darkred',pch=19)
grid()
Figure 6: ex. de la fonction optimize

4.1.2 un exemple de findInterval

rappel de la définition de l’inverse généralisée d’une fonction réelle \(F\), noté \(F^{-1}\)

\[ F^{-1}(x) = \inf\{y\in \mathbb{R}: F(y)\geq x\} \]

x_inv<-findInterval(0.6, cumsum(dbinom(0:10,10,0.2)))
plot(0:10,cumsum(dbinom(0:10,10,0.2)))
points(x_inv,cumsum(dbinom(0:10,10,0.2))[x_inv+1] , col='darkred',pch=19)
grid()
Figure 7: ex. de la fonction findInterval

4.2 TP2

4.2.1 Copule Gaussienne

On commence par un cas particulier (\(dim = 2\)). La copule est un objet pour décrire la dépendance entre les v.a. Pour les v.a. gaussienne, la dépendance est totallement définie par la matrice de covariance.

Soit \((X_1,X_2) \sim \mathcal{N}_2 (0, \Gamma_2)\), où la matrice de covariance est de forme.

\[\Gamma_2 = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\]

Par la construction, on a bien que les lois marginales

\[\mathcal{L}(X_1) = \mathcal{L}(X_2) = \mathcal{N}(0,1)\]

Alors, comme la f.d.r. \(\Phi\in \mathscr{C}^0\)

\[(U_1,U_2) = (\Phi(X_1),\Phi(X_2))\]

est bien la copule qu'on cherche à simuler.

On reviens au cas où \(dim = n\). Si on simule \((Z_1,\dots,Z_n)\) par

\[\forall i \in \{1,2,\dots,n\}\qquad Z_i = \sqrt{\rho} G_0 + \sqrt{1-\rho} G_i\]

où \(G_i\) sont les normals standards indépandantes, c'est facil à vérifier que

\[\mathrm{Var}[Z_i] = 1\qquad et \qquad \mathrm{Cov}[Z_i Z_j] = \rho \qquad \forall i\neq j\in\{1,2,\dots,n\}\]

Alors, la matrice de covariance associée est de la forme

\[ \Gamma_n = \begin{pmatrix} 1 & \rho & \cdots & \rho \\ \rho & 1 & \cdots & \rho \\ \vdots & \vdots & \ddots & \vdots \\ \rho & \rho & \cdots & 1 \end{pmatrix} \]

Donc, pour simuler \((U_1,U_2,\dots,U_n)\) la copule gaussienne associée, il faut juste prendre

\[U_i = \Phi(Z_i) \qquad \forall i \in \{1,2,\dots,n\}\]

L’intérêt de cette méthode est que si vous voulez simuler les vecteurs gaussiens avec \(\Gamma\) différents, on peut simuler \((G_i)_{0\leq i \leq N}\) une seule fois et ceux qui restent sont des calculs numériques. Il nous permet aussi de faire la simulation ‘dépendante’ parallèlement.

4.2.2 Urne d'Ehrenfest

On remarque que la loi invariante est bien \(\mathcal{B}(d,\frac12)\).

\[\forall k \in [0:d], \qquad \mathbf{P}(X_n = k) = \frac{1}{2^d}\binom{d}{k}\]

4.2.2.1 Pourquoi ?

Une manière intuitive de regarder ce problème est de tracer une boule. Imaginons que l’on a déjà fait une infinité de fois de transitions. Alors, une certaine boule soit dans A ou dans B est indépendant que les autres boules, qui suits bien une loi Bernoulli avec proba \(\frac12\). Donc, le cardinal de l’urne A suits bien une loi binomiale.

Or, c’est facile à vérifier que cette chaine est irréductible et récurrente positive (dim. finie). Elle admet une unique loi invariante. Alors il suffit de vérifier que

\[X_0 \sim \mathcal{B}(d,\frac12)\implies X_1 \sim \mathcal{B}(d,\frac12)\]

4.2.2.2 Chaine d'Ehrenfest modifiée

On modifie un peu la modélisation précédente en considérant la règle suivante: “on tire un numéro de balle selon la loi uniforme sur \(\{1,2,\dots,d\}\) et à un tirage \(i\) on déplace la balle numéro \(i\) d’une urne à l’autre avec probabilité 1/2”.

Intuitivement, on a la chaine modifiée est une chaine apériodique irréductible d’états finis, il y a donc une convergence à la vitesse géométrique vers la loi invariante \(\pi\). On note \(\mu_0\) la loi initial et \(M'\) la trasition modifiée

\[\forall f \in \mathcal{B}_b([0:d]),\qquad |\mu_0 (M')^n (f) - \pi (f)|\leq C\alpha^n \Vert f\Vert_{\infty}\]

avec \(C>0\) une const. et \(0<\alpha<1\).

On remarque que la chaine origine est périodique avec période = 2, il n’y a donc pas de convergence vers la loi invariante.

On note \(M\) le noyau de transition pour la chaine origine et \(\pi_0\) la loi invariante associée. on a

\[M'(x, dy) = \frac12 M(x,dy) + \frac12\delta_x(dy)\]

Alors, c’est facil à vérifier que

\[\pi_0 M' = \frac12\pi_0 M + \frac12\pi_0 = \pi_0\]

Par l’unicité de la loi invariante, on a \(\pi = \pi_0\).

4.2.3 Processus de Poisson inhomogène

Soit \(\lambda:[0,+\infty[ \to ]0, +\infty[\) une fonction strictement positive, appelée fonction d’intensité. On considère maintenant des instants de sauts \((T_n)_{n \geq 1}\)

(avec toujours \(T_0 = 0\)) donnés par la dynamique \(T_{n} = T_{n-1} + S_n\) où pour tout \(n \geq 1\), \(S_n\) est une variable aléatoire sur \(]0,+\infty[\) définie conditionnellement à \(T_{n-1}\) par

\[ \begin{aligned} \forall s > 0, \forall t_{n-1} > 0, \quad \mathbf{P} \bigl[S_n > s \; \bigl\vert \; T_{n-1} = t_{n-1} \bigr] &= \exp \bigl(-\int_0^s \lambda(u + t_{n-1}) d u \bigr) \\ &= \exp \bigl(- (\Lambda(t_{n-1} + s) - \Lambda(t_{n-1})) \bigr) \\ \end{aligned} \]

avec \(\Lambda(t) = \int_0^t \lambda(u) d u\) l’intensité intégrée. Dans la suite on note \(\Lambda^{-1}\) la réciproque de \(\Lambda\).

4.2.3.1 Simulation par inverse de l'intensité intégrée

Montrer que si \(E \sim \mathcal{E}(1)\) et que \(t_{n-1}\) est fixé, alors la variable aléatoire \(X = \Lambda^{-1}\bigl(E + \Lambda(t_{n-1})\bigr) - t_{n-1}\) vérifie

\[ \forall s > 0, \quad \mathbf{P} \bigl[X > s \bigr] = \exp \bigl(- (\Lambda(t_{n-1} + s) - \Lambda(t_{n-1})) \bigr), \]

preuve :

Puisque \(\lambda\) est une fonciton strictement positive, l’application \(\Lambda : t \rightarrow \Lambda(t)\) est strictement croissante sur \(\mathbb{R}_+\).

Alors, on a

\[ \begin{aligned} \mathbf{P}\big[ X > s\big] &= \mathbf{P}\big[\Lambda^{-1}\big(E + \Lambda(t_{n-1})\big) - t_{n-1} > s \big]\\ &=\mathbf{P}\big[E + \Lambda(t_{n-1}) > \Lambda(t_{n-1}+s)\big]\\ &=\mathbf{P}\big[E > \Lambda(t_{n-1}+s)-\Lambda(t_{n-1}) \big]\\ &=\exp \bigl(- (\Lambda(t_{n-1} + s) - \Lambda(t_{n-1})) \bigr) \end{aligned} \]

C’est-à-dire, on a

\[\mathcal{L}\big(X \; \big\vert \; T_{n-1} = t_{n-1}\big) = \mathcal{L}\big(S_n \; \big\vert \; T_{n-1} = t_{n-1}\big)\]

□
en déduire une méthode de simulation des instants de sauts \((T_n)_{n \geq 1}\) et du processus de comptage associé (appelé processus de Poisson inhomogène de fonction d’intensité \(\lambda\)).

Tester cette méthode avec la fonction d’intensité \(\lambda(t) = 0.1 + 5t\).

lambda <- function(t) return(0.1 + 5*t)
Lambda <- function(t) return(0.1*t + 2.5*t^2 )
Lambda_inv <- function(s) return(-0.02 + 0.2*sqrt(0.01+10*s))

X_cond <- function(t) {
    E <- rexp(1);
    X <- Lambda_inv (E + Lambda(t)) - t;
    return(X);
}


# svglite("poisson_inhomogene.svg")
Npaths <- 200
T <- 4
plot(NULL, xlim=c(0,T), ylim=c(0,25), xlab="Time t", ylab="Processus de comptage N_t",
     main = paste(Npaths, "trajectoires d'un processus de Poisson inhomogène"))

# on fait simple et on illustre ici l'utilisation d'une boucle for et d'une boucle while
for (k in 1:Npaths) {
    T_n <- 0
    path <- 0
    while (T_n < T) {
        T_n <- T_n + X_cond(T_n) 
        if (T_n < T) path <- c(path, T_n)
    }
    N_T <- length(path)-1
    N <- 0:N_T

    # on ajoute pour faire joli le dernier point (T, N_T)    
    path <- c(path, T)
    N <- c(N, N_T)
    
    lines(path, N, type="s", col=rgb(0,0,1,alpha=0.5))
}
grid()
Figure 8: Processus de Poisson inhomogène

4.2.3.2 Simulation par la méthode des sauts fictifs (thinning)

On suppose maintenant que la fonction d’intensité est bornée

\[ \exists \bar \lambda > 0, \quad \forall t \ge 0, \quad \lambda(t) \le \bar \lambda. \]

Soit \((\bar T_n)_{n \geq 1}\)

les instants de sauts d’un processus de poisson \((\bar N_t) (t \geq 0)\) d’intensité \(\bar \lambda\). On peut simuler les instants de sauts \((T_n) (n \ge 1)\) en supprimant certains sauts de \((\bar T_n) (n \ge 1)\). Plus précisément, on considère \((U_n)(n \ge 1)\) une suite i.i.d. de loi uniforme sur \([0,1]\) et on construit la suite d’indices \((\tau_n)(n \ge 1)\) par la récurrence

\[ \forall n \ge 0, \quad \tau_{n+1} = \min \Bigl\{ k > \tau_n, U_k \le \frac{\lambda(\bar T_k)}{\bar \lambda} \Bigr\} \qquad \tau_0 = 0 \]

Alors les sauts de \((\bar N_t)(t \ge 0)\) d’indices \((\tau_n)(n \ge 1)\) sont les sauts d’un processus de Poisson inhomogène de fonction d’intensité \(\lambda\). On peut alors construire \((T_n)(n \ge 1)\) en posant \(T_n = \bar T_{\tau_n}\).

tau_update <- function(tau,lambda_bar, T_bar){
    
    # tau: tau_n
    # lambda_bar: la borne de l'indensité
    # T_bar: un vecteur des instants de sauts d'un processus de Poisson
    
    k <- tau
    while(runif(1) > T_bar[k]/lambda_bar){
        k <- k+1
    }
    # la sortie: tau_{n+1}
    return(k)
}

4.2.3.3 Variable de contrôle

Soit \(X\) une v.a. et on s’intéresse à estimer sa moyenne \(\mathbf{E}[X]\) par la méthode de monte-carlo. Soit \(Y\) une v.a. telle que l’on sait la moyenne \(\mathbf{E}[Y]\). Alors, on a

\[\mathbf{E}[X] = \mathbf{E}[X-\alpha Y] + \alpha \mathbf{E}[Y]\]

On suppose que on sait simuler suivant la loi \(\mathcal{L}(X)\). Alors, on s’intéresse au cas où

\[\mathrm{Var}[X - \alpha Y] \leq \mathrm{Var}[X]\]

C’est facil à vérifier que une condition nécéssaire est que

\[\mathrm{Cov}[X,Y] \neq 0\]

Et si \(\alpha \in (0,\frac{-2\mathrm{Cov}[X,Y]}{\mathrm{Var}[Y]})\) (ou \(\alpha \in (\frac{-2\mathrm{Cov}[X,Y]}{\mathrm{Var}[Y]},0)\)), on a \(\mathrm{Var}[X - \alpha Y] \leq \mathrm{Var}[X]\).

De plus, \(\mathrm{Var}[X - \alpha Y]\) atteint son minimum lorsque \(\alpha = -\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}[Y]}\).

4.3 TP5

4.3.1 Algorithme de Métropolis-Hastings

Ici, on donne une version un peu plus générale que dans la feuille de TP : On utilise opérateur d’intégrale au lieu de matrice de transitions. Le but est de comprendre pourquoi M-H nous donne un noyau de transition avec la bonne loi invariante et comment on fait la simulation. On remarque que tous les calculs seront identiques si on utilise la représentation matricielle. On va d’abord voir la construction de noyau de transition et vérifier qu’il nous donne la bonne loi invariante par les calculs algébrique. En suite, on donne deux décompositions pour mieux comprendre le dynamique associé et l’idée de la construction. On remarque que ce genre de dynamique nous permet de construire les autres algo pour les sénarios différents.

On suppose qu’on a une mesure \(\pi(dx)\) qui est absolument continue par rapport à mesure \(dx\), disons mesure de Lebesgue et un noyau de transition \(M(x,dy)\) tel que pour tous les \(x\), \(M(x,dy) = M(x,y)dy\) est une mesure de proba et est aussi absolument continu par rapport à \(dy\). On suppose que \(M\) est \(\pi\)-reversible, c-à-d

\[ \pi(dx)M(x,dy) = \pi(dy)M(y, dx) \] On remarque que l’égalité est au sens que pour tous les fonction de test \(f\) bornée mesurable, on a

\[ \int f(x,y) \pi(dx)M(x,dy) = \int f(x,y) \pi(dy)M(y, dx) \] Alors, on définit le noyau Métropolis-Hastings par

\[ Q(x,dy) := r(x,y)M(x,dy) + \left(1 - \int r(x,s)M(x,ds) \right)\delta_x(dy) \] avec \[ r(x,y) = \frac{\pi(y)M(y,x)}{\pi(x)M(x,y)} \]

C’est facil à vérifier que \[ \int f(y) \pi(dx)Q(x,dy) = \int f(y) \pi(dy) \] pour n’importe quelle fonction de test mesurable bornée \(f\).

4.3.1.1 Décomposition 1

Vu que \(r(x,y)M(x,dy)\) n’est pas un noyau markovien (\(\forall x, \int M(x,dy) = 1\)). Donc, on en déduit que

\[ Q(x,dy) := \int r(x,s)M(x,ds) \cdot \frac{r(x,y)M(x,dy)}{\int r(x,s)M(x,ds)} + \left(1 - \int r(x,s)M(x,ds) \right)\delta_x(dy) \] où \[ \frac{r(x,y)M(x,dy)}{\int r(x,s)M(x,ds)} \] est bien un noyau markovien. Alors, pour chaque état \(x\), on pourrais regarder la loi de \(y\) comme une loi de mélange de \(\frac{r(x,y)M(x,dy)}{\int r(x,y)M(x,dy)}\) et \(\delta_x (dy)\), avec proba \[ \left( \int r(x,s)M(x,ds) , 1- \int r(x,s)M(x,ds)\right) \]

4.3.1.2 Décomposition 2 (pour simulation)

En même temps, on a une autre décomposition pour simplifier la procédure de simulation. On considère deux étapes :

  • simuler \(\widetilde{Y} \sim M(x, ds)\)
  • simuler \(u \sim \mathscr{U}(0,1)\)
    • si \(u < r(x, \widetilde{Y})\), on pose \(Y = \widetilde{Y}\)
    • sinon, on pose \(Y = x\)

C’est équivaut à dire que on prend

\[ \widetilde{Q}(x,dy) := \int_s M(x, ds)\left\{ r(x,s)\delta_s(dy) + (1-r(x,s))\delta_x(dy) \right\} \] C’est aussi très facil à vérifier que \[ Q(x,dy) = \widetilde{Q}(x,dy) \]